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Abstract. Axion dark matter can resonantly convert to photons in the magnetosphere of 
neutron stars, possibly giving rise to radio signals observable on Earth. This method for 
the indirect detection of axion dark matter has recently received significant attention in the 
literature. The calculation of the radio signal is complicated by a number of effects; most 
importantly, the gravitational infall of the axions onto the neutron star accelerates them 
to semi-relativistic speed, and the neutron star magnetosphere is highly anisotropic. Both 
of these factors complicate the calculation of the conversion of axions to photons. In this 
work, we present the first fully three-dimensional calculation of the axion-photon conversion 
in highly magnetised anisotropic media. Depending on the axion trajectory, this calculation 
leads to orders-of-magnitude differences in the conversion compared to the simplified one- 
dimensional calculation used so far in the literature, altering the directionality of the produced 
photons. Our results will have important implications for the radio signal one would observe 
in a telescope. 
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1 Introduction 


Axions comprise a class of hypothetical particles that were first conceived as a consequence 
of the ‘Peccei-Quinn’ solution to the strong CP-problem in particle physics [1-4]. Axions are 
naturally produced in the early universe, and provide an increasingly popular candidate for 
explaining dark matter [5-12].! Testing the axion dark matter hypothesis is a central goal 
of contemporary fundamental physics. 

Any search for axion dark matter relies on the transfer of energy from the astrophysical 
axions into Standard Model particles. One promising route to achieve this is through the 
axion-photon coupling, which receives ‘universal’ contributions from axion-pion mixing, thus 
providing a rather robust and model-independent target. Most experimental and observa- 
tional ideas to search for axion dark matter are based on the axion-photon coupling. 

Astronomical searches for axion dark matter via the axion-photon coupling seek to 
leverage strong astrophysical magnetic fields to effectively transfer energy from the axion 
field into photons (or vice versa). Particularly interesting are environments where plasma 
effects generate an effective mass for the photon that allows the axion and photon dispersion 
relations to intersect. In such a degenerate situation, dark matter axions can resonantly 
convert into photons, enhancing the photon signal and possibly leading to signals detectable 
with telescopes on Earth. 


'For the purposes of this work, we indiscriminately refer to both QCD axions and axion-like particles 
(ALPs) when we write ‘axion’. For our numerical results, we will assume a relation between the axion mass 
and the axion-photon coupling as appropriate for a QCD axion, but in all of our equations we keep the axion 
mass and its couplings as independent parameters allowing to apply our results to both QCD axions and 
ALPs. 


Neutron stars are ideal astrophysical environments for converting dark matter axions 
into photons [13, 14]: Neutron stars can host extremely strong magnetic fields (up to ~ 
10'° G), and are surrounded by a plasma that admits resonant axion-photon conversion near 
the surface of the star. Due to the plasma densities in the neutron star magnetospheres, 
resonant conversion is most promising for axions in the neV—meV mass range, including the 
perhaps best motivated w~eV—meV mass range for QCD axion dark matter. The corresponding 
photon signal would be a narrow spectral line in the MHz-THz (radio) frequency range, 
with its frequency set by the axion mass. Radio telescopes combine excellent spatial and 
frequency resolution with sensitivity to very weak radio signals. For these reasons, axion- 
photon conversion in neutron star magnetospheres has recently become a very active topic of 
research, and significant effort has been devoted to characterising the possible signal [15-24], 
and searching for it in radio data [25-28]. 

In order to calculate the radio signal from the conversion of dark matter axions into 
photons in the magnetosphere of neutron stars, one must model a number of effects. First, the 
phase space of the axion dark matter evolves under the gravitational influence of the neutron 
star. Second, in the region where the effective photon mass approximately matches the axion 
mass, one needs to solve the system of equations describing the axion-electrodynamics to 
compute axion-photon conversion. Finally, one must propagate the electromagnetic modes 
excited by the axion through the inhomogeneous magnetosphere to the observer far away from 
the neutron star. Regarding the first effect, previous work has either assumed that axions are 
on purely radial orbits (in the frame on the neutron star) as they cross the zone of resonant 
conversion, or that the axion velocity distribution is isotropic at the conversion zone; neither 
of these assumptions is correct. We discuss how to compute the appropriate axion phase- 
space density under the influence of the gravitational field of the neutron star for a neutron 
star with speed relative to the rest frame of the dark matter distribution comparable to the 
dark matter velocity dispersion, as one expects for a realistic neutron star. The main work in 
this paper is to present a calculation of the resonant axion-photon conversion in the neutron 
star magnetosphere (i.e., the second effect above) which fully accounts for the anisotropic 
plasma density in the magnetosphere. This fully three-dimensional calculation improves 
over the simplified one-dimensional calculation previously used in the literature, and we will 
discuss it in more detail in the following paragraph. We do not perform any calculation of 
the third effect, the propagation of the photons from the conversion zone to the observer at 
infinity, in this work. This effect can be dealt with by ‘ray-tracing’ calculations, following the 
electromagnetic modes through the anisotropic plasma of the neutron star magnetosphere, 
see Refs. [20, 23, 24]. 

Let us now discuss the core of this work, the calculation of resonant axion-photon con- 
version in highly magnetised anisotropic plasmas, e.g., in the magnetosphere of a neutron 
star. A central assumption of essentially all previous studies of this problem is the use of an 
effectively one-dimensional formalism for calculating the conversion rate of axions into pho- 
tons, based on the well known calculation for transverse modes [29].? The one-dimensional, 
linearised system can be expressed as a ‘Schrédinger-like’ equation, with time replaced by 
the spatial coordinate along the direction of propagation. The ‘conversion probability’ of ax- 
ions into photons calculated from this equation can then be interpreted as the flux transfer: 


? An interesting partial exception is Ref. [19], which acknowledged the limitations of the one-dimensional 
formalism and developed several results for a two-dimensional setting, but left the full three-dimensional 
problem for further studies. Similarly, in Ref. [23] the change of media in other directions was taken as a 
limiting factor in the conversion probability without solving the 3D equations. 


i.e. the ratio of energy densities of the transverse mode of the electric field and the axion. 
However, in this paper we show that the intuition built around transverse modes breaks down 
in highly aniosotropic media, such as the magnetosphere of a neutron star. The anisotropy 
causes a strong mixing between different photon polarisations, so that propagating states are 
no longer purely transverse or longitudinal.® It follows from this that the naive ‘conversion 
probability’ no longer describes the physically interesting flux transfer. 

We revisit the classical problem of axion-photon mixing in anisotropic media, including 
relativistic plasma effects, and derive a new ‘Schrodinger-like’ equation governing the evo- 
lution of the Langmuir-O (LO) mode, which is the propagating mode in the plasma that 
is excited by axion-photon conversion. Our calculation has several important implications. 
First, the anisotropy of the neutron star magnetosphere results in the photon mode evolving 
along a different direction than that given by its momentum. This invalidates the one- 
dimensional analysis used in previous work, and can lead to orders-of-magnitude change in 
the axion-photon mixing for certain axion trajectories. Second, the ‘Schrédinger-like equa- 
tion’ now includes damping/growth terms, that can invalidate the unitary evolution of the 
system. Third, the energy stored in the LO modes is not given by the transverse portion of 
the field only, and the full Hamiltonian must be used to calculate the flux transfer. Con- 
sequently, the flux transfer is no longer identical to the ‘conversion probability’ calculated 
from the Schrodinger-like equation, so that a naive over-reliance on the quantum mechanical 
analogue would lead to incorrect conclusions about axion-photon mixing. 

The remainder of this paper is organised as follows. We begin by discussing axion- 
photon conversion in an anistropic plasma and a magnetic field in section 2, and derive 
the new form of the ‘Schrédinger-like’ equation. In section 3, we find the relation between 
the ‘axion-photon conversion probability’ and the flux transfer. In section 4 we discuss the 
calculation of the axion phase space in the neutron star magnetosphere for an neutron star 
moving with arbitrary speed relative to the rest frame of the dark matter distribution. In 
section 5 we present numerical results for our conversion calculation and compare them to 
the results one would obtain in the simpler one-dimensional calculation. As we will see, the 
conversion probability can differ by orders of magnitude between our calculation that fully 
accounts for the anisotropy of the plasma and the simplified one-dimensional calculation. We 
conclude in section 6. Further technical details of our calculation of axion-photon conversion 
can be found in the appendices. 


2 Axion-photon conversion in aniostropic plasmas 


We consider dark matter axions that are gravitationally accelerated towards a neutron star 
and may convert to photons in its magnetosphere. Due to the large gravitational field of 
the neutron star, the axions will generally be mildly relativistic, a particularly complicated 
regime from the point of view of calculating the axion-photon conversion. In such a regime, 
medium effects cannot be ignored, however the axion also has significant momentum and so 
the conversion has a strong directional dependence. Our calculation will focus on the regime 
where the WKB approximation can be applied, i.e. where the axion has momentum larger 
than the first derivatives of the E-fields. As long as this condition holds, our calculation 
will apply regardless of the relativistic (or not) nature of the axion. Axion-photon mixing is 
maximised in regions satisfying the resonance condition, which is approximately given when 


3While this has been noted in [23], the implications for how one calculates axion-photon conversion have 
not yet been explored. 


the plasma frequency is equal to the axion mass wp ~ Ma. Since this condition is typically 
fulfilled only in small regions, we solve the mixing problem in the flat-space approximation, 
neglecting gravity. This approximation is expected to be accurate for most (but not all) 
possible trajectories of infalling axions. 

The classical theory of axions coupled to electromagnetism is governed by the Klein- 
Gordon equation and Maxwell’s equations, modified to include the axion-photon mixing. 
Following common practice (see, for example, Refs. [14, 30, 31]) we can linearise the equa- 
tions of motion in a background with the magnetic field Byg, and restrict (without loss of 
generality) to electromagnetic modes that oscillate with frequency w: 


wat V2a = ma — gayE- Bys (2.1) 
-V?E+V(V-E) = w°D + w’ gaya Bys, (2.2) 


where we have set the relative permeability to u = 1. The electric field is denoted by E, and 
the displacement field by D. We would now like to solve these equations for dark matter 
axions traversing a neutron star magnetosphere. 

We must first express the equations (2.1) and (2.2) in a convenient coordinate system 
for any given axion. Two basis vectors can be defined respectively by the direction of motion, 
and the external magnetic field, which generically are linearly independent. We denote the 
axion’s direction of motion by Z, so that k = kz. This will also be the initial direction of 
any produced photons, though the changing medium will result in the photon momentum 
evolving over longer distance scales. We take Byg to lie in the (y, z)-plane at an angle 0(z) 
with the positive z-axis. In the following, we assume 0 4 {0,7} unless otherwise specified. 
This defines (z) and consequently also X(z) by X(z) = (z) x 2. We will see that the 
conversion generally occurs over a small region, and so we will neglect gravity and assume 
that the axion travels along a straight line and £4 = 0. It then follows that 4g x y and 
Ly œ &. The expressions for the differential operators V? and V(V -E) in the (£, 7,2) basis 
are then identical to those in a fixed Cartesian coordinate system. In general, 0 varies due 
to gravitational bending, the resulting photon refracting in the medium and the fact that 
the magnetic field itself changes. Again as the region where axion-photon conversion occurs 
should be small, we take 0 to be almost constant. 

The electric displacement field for a highly magnetised plasma is given by the constitu- 
tive relation [32] 

1 0 0 
D= ej EI —{01- (ser) sin? 6 (3 ) cos 8 sin 0 -| E 


G2 
2 2 
0 (ser ) cos Osin 6 1 — (38) cos? 6 Ez 


(2.3) 


where ¢€;; is the dielectric tensor with the plasma frequency wp. One can see that the 
anisotropy of the dielectric tensor mixes E, and E,, so one is not in general left with purely 
transverse of longitudinal propagating states. Here, we have assumed that the particles form- 
ing the plasma are moving along the field lines with speed 1) in the neutron star frame, and 
define i 

=w- kup y= ——— =. (2.4) 
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We have also defined T 
(..) = ay dp) ---F}\ 1) (2.5) 


where Fi is the distribution function of the species 7 in the plasma. Note that for an ultra- 


relativistic plasma, (773) ~ (y)~! [33]. In Refs. [14, 19, 20], the electrons in the plasma were 
treated as non-relativistic, i.e. y = 1. While we allow for relativistic electrons, we note that 
infalling axions will be at most mildly relativistic. For mildly relativistic axions, one would 
expect © ~ w, however, in principle the distinction between w and w can lead to differences 
for high momenta axions traveling through relativistic plasmas. For dark matter falling into 
a neutron star, this will simply change the effective wp by a factor of @/w which is trivial 
compared to the significant modelling uncertainties of neutron stars (discussed in section 
oe 

When the axions are decoupled from the system, one can solve for the propagating 
electromagnetic states of the medium. Note that such states only apply locally: as one 
moves through the medium, the eigenstates deform adiabatically. As noted in [23], the 
relevant electromagnetic mode (i.e. the mode excited by axion-photon conversions) is the 
so-called Langmuir-O (LO) mode, whose polarization has both longitudinal and transverse 
components. As we discuss in Appendix A, in principle, the axion can also excite purely 
longitudinal modes. However, such purely longitudinal modes are non-propagating and are 
hence irrelevant for the electromagnetic signal observed far from the neutron star. The LO 
mode, on the other hand, evolves adiabatically into a purely transverse vacuum mode as 
it exits the neutron star atmosphere [33]. The specific form of the LO mode is found in 


Appendix A to be 
z —1 : a cos @ sin 0 : 
Ero = y 2 , (2.6) 


= Z 
1 D4 cos? 0 sin? 0 wt — w2 cos? 0 
+ (27—02 cos? 0)? 


where we have defined an effective relativistic plasma frequency @? = ww? /@*). In the 
non-relativistic limit, as has been mostly considered in the literature, there is no distinction 
between Wp, and wp. Our task is to calculate the conversion of axions into the LO mode. While 
the modification to the dispersion relation was noted in Ref. [23], the effects of the mixing 
of transverse and longitudinal modes on the conversion calculation and the implications for 
how one defines a photon were not explored. As we will see, such effects can change the 
conversion probability for a given axion trajectory by orders of magnitude. 

We can now write out the z, y, z components of equation (2.2). We will be considering 
the case of axion photon mixing in a slowly varying plasma, so we can neglect second order 
derivatives that do not involve the z-direction. In addition, EF, does not directly couple to 
the axion as, by construction ¢1Byg. While Ey can be generated by derivative terms, it is 
suppressed when compared to Ey > so we can also neglect it. At lowest order this leaves the 
following coupled equations: 


PE, 2 ee 2 OPE. 2 

Deby x (w — wp sinf 0) Ey + ©% cos 0 sin 0E; + aE — w Jaya Bys sin 6 , (2.7) 
2 

oy ~ (w? — ©? cos” 0) E; + ©? cos 0 sin 0Ey + w’ gaya Bys cos 6. (2.8) 

OzOy L i 


Previous studies of axion-photon conversion near neutron stars have all neglected the left- 
hand side of these equations, featuring mixed y- and z-derivatives.* Neglecting mixed deriva- 


“While it is noted that the full gradient needs to be dealt with in Ref. [19, 23], the conversion was not 
calculated including mixing terms. All previous calculations including Ref. [14, 20] similarly neglect the mixed 
derivatives. 


tives ignores the mixing between E, and E, and is equivalent to assuming that Ey only 
evolves along z. This does not hold for general trajectories, and is only a good approxi- 
mation if 0 ~ 7/2, i.e., when the axion happens to be traveling perpendicularly to Bys. 
To simplify matters, we can eliminate E, from the equations. This will allow us to solve 
directly for the propagating mode, and neglect non-propagating contributions to E,, which 
are discussed in Appendix A. By substituting equation (2.8) into (2.7) and neglecting any 
higher order derivatives that follow, we end up with a single partial differential equation for 
Ey, and a: 


OE, 2&2 cosdsind GE, 2w*dpysinOcosO AH, OEy  2w*GpsinOcosO OB, OEy 
O?z w? —G2cos?@0zdy (w? — w? cos? 0)? Oy Oz (w? — w2 cos? 4)? Az dy 
w — a? w? sind 
~ P? E Jaya Bys . (2.9) 


02 y @2 
1— 4 cos? 0 I= cos? 0 


We neglect also the derivatives of the axion driving term, as they will only have a small 
impact on the evolution of Ey. 

To attempt an analytic solution, we assume that the change in @,, and so k, is small 
and thus the direction of propagation stays the same. Of course, this is not true in gen- 
eral: photon trajectories will bend in a changing medium. Assuming small changes in W,, k 
necessarily restricts us to solving over some finite region where the photon path is approxi- 
mately constant. The region over which the majority of the generation of an axion-induced 
E-field occurs is referred to as the conversion zone. While one can estimate the radius of 
curvature for the changing direction of the photon, such as was done in Ref. [23], our focus 
in this work is not on the propagation of the converted photons. As most of the generation 
of the propagating LO mode via the axion should be localised to a smaller region, separat- 
ing propagation and conversion allows for an analytic approximation of the fields generated 
by the axion, which could then be propagated numerically via ray tracing code such as in 
Refs. [23, 24]. 

As long as we are treating the medium as slowly varying, we can write the fields as an 
envelope term, denoted by a tilde, multiplying a plane wave, 


alejet» | (2.10) 


Ey = Ey(y, geek) , a 


To move further, we will take the WKB approximation and assume that kĒy(y,z) > 
OE, (y, z)/Oy, OB, (y, z)/dz > OE, (y,z)/dz?. As we will see, the smaller the derivatives 
the longer the length scales over which axions can resonantly convert, leading to higher con- 
version probabilities. Thus a breakdown of the WKB approximation, which occurs when 
the first derivatives of the @, vanish, corresponds to a region of enhanced axion-photon 
conversion. The WKB approximation allows us to then simplify Eq. (2.9) to 
a? J 5 2 
5 A a ~ (m? — £02 — ikD) By - amA (2.11) 


Qik 
R sin 


where we have used the on-shell relation for the axion, m2 = w? — k?, and defined 


sin? 0 


— on” 
(73) 
1 — Z cos? 0 
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(2.12) 


The damping/growth term D appears as an artifact of assuming that the photon does not 
experience any curvature, i.e., that it is well described as only having momentum in the 
z-direction, and is given by 
Qup£? Oy 
wtan@sin? 0 Oy ` 
The back-reaction on to the axion field should be small, so we can treat à as constant. 
To progress, we can write a new differential operator 
-2 

ane + we 0 (2.14) 

ðs Oz w*tan@ dy 
This operator allows us to finally write the axion-photon mixing equations as a Schrödinger- 
like equation 


(2.13) 


(m? — £02 — ikD) Ey - 


OB, 1 1 we 

joe = ăBys, 2.15 
Os 2k Dk sin 9907 NS ( ) 

As 0/0s is a function of ùp and 0, and in general &p(s), any integration over s will in principle 

be curvilinear. However, if we are only concerned with small regions, we can assume that œp 

is only slowly varying and so s is constant with respect to y, z. Assuming a constant W, we 


can write =? 
1 cos 0 sin 0 
S~ a Peace y+a). (2.16) 
1 ws cos? 0 sin? 0 Wp COS 0 
+ (w?—W2 cos? 0)? 
Comparing with Eq. (2.6) we see that § is actually perpendicular to Exo. In other words, 
the propagating mode Fro is transverse with respect to the direction it evolves in, rather 
than the direction of propagation. While the mode is not transverse, it evolves in a similar 
manner to a transverse wave. Note that 0/0s is not defined with unit norm for simplicity, 


however that means one has to be slightly careful when talking about length scales associated 
with s. For non-relativistic axions (Wp = Ma), we can write down simple expressions 


Bug = —sin0¥ + cos0z, (2.17a) 
Exo ~ -sinb ŷ + cos 04, (2.17b) 
S§~cosd¥+sindzZ. (2.17¢) 


In this limit we can simply show the relationship between Bys, z, ELo and s, depicted in 
Fig. 1. We can see in this limit that Bys||ELo: the LO mode couples perfectly to the axion 
via E. B. For ultrarelativsitic axions with w >> @,, we see instead Ê ~ —ŷ, 8 ~ ĉ and 
E œ~ sin? @ giving back the usual mixing of an axion with a transverse photon [29]. 

Ignoring the overall phase, as shown in Appendix B, Eq. (2.15) can be solved in the 
stationary phase approximation to give 


E,|so + L/2] = w? il elo” ds'D/24, Bys , (2.18) 
2k [epas + w? L 36) 


where &, = ĝwp/3s, we have allowed for a small change in 0 given by 6’ = 00/0s and so is 


given by solving for the resonance condition 
meu 


Ge) = (2.19) 


m2 cos? 6 + w? sin? 0 
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Figure 1. We show the various axes associated with axion-photon conversion in an anisotropic 
medium. The axion propagates in the Ê direction, which is at an angle 0 from the magnetic field 
Byg. The magnetic field is defined to lie in the z,y plane. In the non- relativistic limit, the E-field 
associated with the propagating LO mode, Exo, is given by Bys, ice. , Exo is parallel to the B-field. 
The magnitude of Exo evolves along the perpendicular direction, denoted by ŝ. 


One can see how the LO mode extrapolates between longitudinal and transverse behaviour 
as one varies 0. If 0 = 7/2, the dispersion relation is that of a regular transverse mode, with 
a resonance given by Ùp = Ma. If instead 0 = 0, then the resonance condition is that for 
an axion converting to a longitudinal plasmon with negligible spatial dispersion, ©, = w, as 
for example studied in Ref. [34]. In the non-relativistic limit, w — Ma, and the resonance 
condition simply becomes Wp» = Ma for any 0. The conversion length can be read off from 
the exponent in the stationary phase approximation (similar to Ref. [19]), giving 


sin 0 tk tk 


L= we a ~ sinl Za (2.20) 
£ Opi, + roo" Wp|a, 


where the latter statement holds in the non-relativistic limit. To evaluate the damping term 
in Eq. (2.18), we note that by assumption, the equation is only valid in a small region 
around the conversion zone where higher order derivatives and the bending of photons can 
be neglected. Assuming that 0W,/Oy = O(0W,/0s), one can see that over the conversion 


zone (so E L/2) 

so+L/2 k 

f ds'D = O(a 5 ) <1, (2.21) 
so—L/2 aL 


thus, for simplicity, we will neglect this damping in the probability of conversion. We can 
also include the longitudinal component of the LO mode, giving 


P p T 

Eylso + L/2] ~ gayaBngw* T (2.22) 
2k (opa, o =e E 20 
g2 

7 cos 0 sin 0 

BE, [59 + L/2] ~ gayiBygw?—?—— — l (2.23) 
w2 — w2 cos? 0 wp 

P 2k Wp + 20 


Actually, as discussed in Appendix A, there are two contributions to Æ. As we are concerned 
with propagating modes that may exit the neutron star, here we only consider the longitudinal 
component of the LO mode. There is also a term that directly comes from the axion-E, 


mixing, however, it only persists in the region where the axion mixing is strong, and follows 
the phase of the axion, rather than the LO mode. Thus, after a length L, this axion-mixing 
term will dephase with the LO mode and not contribute to the E-fields exiting the neutron 
star. If we then take ùp = w = mg, we get the simple expressions 


3 
~ . | nm 
Ey [so +L/2] = Jay& Bys AA] j (2.24) 


? m2 cos@sin@ | mm 
E.[s9 + L/2] © gayaB = 7° 
z[so + L/2] = gay4Bns k2 + m2 sin? 0 \/ 2k|a%,| 


(2.25) 
We kept the factor of k? in the denominator of E, to avoid an unphysical divergence at 
6 — {0,7}. Note that the E, component can actually dominate the E-field: 


E, mēcosðsinð 
Ey  k2+m?2sin?6’ 


(2.26) 


which has a maximum at |E;/Ey| ~ mq/2k for k < ma. Even for large k, Ey and E; are 
generically of similar magnitude. 


3 Probability of conversion 


When expressing the classical problem of axion-photon mixing as a Schrédinger-like equation 
(with time replaced by the spatial coordinate z or s), it is natural to make use of the language 
and techniques originally developed within quantum mechanics to simplify calculations. The 
most natural object to compute from the Schrödinger equation is the ‘conversion probability’ 
(|Ay/a|?), i.e. the norm square of the transition amplitude of going from a pure axion state 
(a = 1) into a pure photon state (A,). However, it is important to remember that this is 
not actually a quantum mechanical calculation, for which one would have to use the correct 
Schrödinger equation (with time derivatives and h 4 0), and carefully defined wavefunctions. 
In the vaccuum case, the ‘conversion probability’ calculated from the classical Schrödinger- 
like equation is of direct physical interest as it gives the ratio of the energy densities stored 
in axion and electromagnetic fields [29]; this definition of conversion probability has been 
used in Refs. [14, 19, 20, 23] in the context of axion-photon conversion near neutron stars. 
However, as we discuss in this section, when taking into account the 3D complexity of the 
neutron star magnetosphere and the correct definition of propagating modes in the plasma, 
the naive ‘conversion probability’ is no longer the physically relevant object to compute. In 
essence, such a definition does not include the fact that the Hamiltonian is modified by the 
presence of matter, so would miscount the number of photons generated if interpreted as an 
actual conversion probability. 

The key point of this section is that the (time-averaged) energy density stored in the 
propagating LO mode is no longer simply 1/4(1 +o /w?)|E,|?, as is the case for an isotropic 
non-relativistic plasma. To calculate the energy stored in the LO mode, we must consider 
that the medium has both temporal dispersion, and is anisotropic, giving [35, 36] 
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1 
= 480! 1 


wed) E} E; + i 


HZH’ . (3.1) 


While the H-field terms are suppressed at order k?/w? relative to the E-field terms, near 
the neutron star surface, k can be relatively large compared to w, meaning that significant 
energy can be stored in the H-fields. For a propagating state given by E,o, we find 


(Qu* — wc? + (o5 — 3w 0) cos(20) + ©) 


U = |Ey|? 2 
8 (w? = w2 cos? 0) 


2 k? 
+ [Ey (3.2) 


Equation (3.2) reduces to 1/4(1+@?/w?)|E,|? in the non-relativistic limit and when 8 = 7/2, 
and reduces further to |Ey|?/2 for w = Wp, and it is only in this case that the naive conversion 
probability (|A,/a|”) gives the ratio of energy densities stored in the electromagnetic and 
axion fields. Notably, the energy stored in the propagating mode is enhanced for B-fields 
almost aligned with the direction of propagation. In the absence of medium losses, the energy 
stored in the photon field should be conserved: when propagated through a slowly changing 
medium the propagation state will adiabatically deform, but maintain stored energy. Thus, 
taking only the transverse component of the stored energy into account, as done in previous 
studies, neglects a significant part of the stored energy, or in other words, miscounts the 
number of photons that are converted. 

For a simple expression, we can finally write the ratio of the energy densities of the 
electromagnetic propagating mode and the axion field in the non-relativistic limit (ma > k) 
as 


a 2U Ae (2wt w? a | (cs 3w* a?) cos(20) + ws) P k2 
wal? a? 4 (w? = we cos? 9)” 2w? 
2 p2 
B 5 
w Jar ?Ns uk sin? 8. (3.3) 


a 
~ 2k|õl| (k? + m2 sin? 6)? 


As the system is considered here to be time independent, a single axion should convert 
to a single photon with energy w and so the ratio of energy densities should correspond 
to the ratio between the axion and photon fluxes. Note that this flux transfer is distinct 
from the ‘conversion probability’ (|A,/a|*) calculated from the Schrédinger-like Eq. (2.15), 
essentially reweighting it via the Hamilitonian. However, if one did a full quantum mechanical 
calculation, the expectation value of the conversion probability of axions to LO photons 
should agree with this flux transfer [37]. 

Unlike the result of Ref. [14], no propagating photon is generated from a non-relativistic 
axion in a longitudinal B-field (i.e., for Byg||Z). However, we should stress that 0 = {0,7} 
is a special condition and the order of limits matters. For purely longitudinal conversion one 
should perform a dedicated calculation (i.e., solve Eq. (2.8) directly or perform an analysis 
in the vein of Ref. [34]). We should note that Eq. (3.3) breaks down under more extreme 
conditions. If |@,| = 0, in other words, if @p is constant along s, the conversion length 
becomes infinite. To regulate such divergences, one must cut-off L when it is competitive 
compared to other limiting scales, such as the radius of curvature of the photon trajectory 
or that over which second derivatives are relevant [23], cf. Appendix B. 

Our ratio can be contrasted with the ‘conversion probability’ calculated when the non- 
trivial Hamiltonian and the mixing of transverse and longitudinal modes are neglected, 


pid TW 2 


= —— __? Beasin 0. A 
ay 2k|ðp/3z| I Ns SIN (3 ) 
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If we take 6 = 1/2 so that 8 = z, reducing back to a purely one dimensional problem, then 
we see that both expressions agree, as one would expect 


T TW 
R> PIR, (0=5) = G2 Bis. 3.5 
ay 2 kja 27 NS ( ) 
We have now derived in Eq. (3.3) an analytic expression for axion-photon conversion in 
highly aniosotropic systems. In cases where the axion is traveling at a non-trivial angle to 
the external B-field there can be significant alterations to the flux transfer. To see what 
impact the 3D calculation has in a neutron star, we must consider a specific neutron star 
model. 


4 Phase space distribution 


In this section, we discuss the phase space distribution of the axions at the conversion surface. 
An excellent discussion of the change of the phase space distribution of collisionless particles 
under the influence of the gravitational field of an astrophysical object such as a neutron 
star can be found in Ref. [38]. We will assume that far away from the neutron star, the 
velocity distribution of the axions in the galactic rest frame is given by a Maxwell-Boltzmann 
distribution with velocity dispersion o, truncated at the Galactic escape velocity Vege as in 
the Standard Halo Model [39-41], 


Pies ( l J EEE (oe =o) , (4.1) 


Nese \ 2702 


where H (x) is the Heaviside step function, and the normalization factor is given by 


2 
Vesc 2 Vesc — fee 
Nesc = E 4/ 205 , 4.2 
esc rr ($=) T Oy e ( ) 


We can shift this distribution into the frame of the neutron star, 


fv) = R (v — vns) , (4.3) 


where vyg is the velocity of the neutron star relative to the galactic rest frame. The six- 
dimensional phase space density far away from the neutron star is then 


IPO) = n& for(v) = PE for(w) , (4.4) 


where n2° (pẹ) are the axion number (mass) density far away from the neutron star. 
Liouville’s theorem states that the phase space density is conserved along trajectories. 
Hence, at some point rng in spherical coordinates centered on the neutron star, the phase 
space density is given by 
felens; v) = [2° (væ) ; (4.5) 


where Voo = Vo(ryns, V) is the velocity at infinity for an orbit with velocity v at ryg. To 
distinguish from our axion-centric frame as used in the previous sections, we use NS as a 
subscript; in other words, the position vector is given by fns = (rys, Ong, Ons) in spherical 
coordinates. Let us stress here that Liouville’s theorem states only that the phase space 
density is conserved along trajectories, it does not make any immediate statements about 


= Jl = 


the relation of the functional form of the phase space density at different r. In order to 
obtain that functional form, one must compute voo(rns,v). For a Newtonian potential, 
energy, angular momentum, and the Laplace-Runge-Lenz vector of a particle moving in the 
potential are conserved, such that [38] 


v2 v + væ (GMns/TNs) ENS — vov (v - rys) 
v2, + (GMys/ryns) — vo (V - rys) 


(4.6) 


VælrNs, V) = 


with vœ = J v? — 2GMnzg/rTng, the gravitational constant G, the mass of the neutron star 
Mys, and fns = rns/TNs- 

It is interesting to note that due to the symmetry of the system, if the velocity and 
mass distribution of the axions is homogeneous and isotropic (in the frame of the neutron 
star) far away from the neutron star, i.e. vns = 0 in Eq. (4.3), the velocity distribution will 
be isotropic at any rys. In this case, one can immediately compute [20] 


foltns,¥) = © > ( 


Ma Nese 


1 
2no? 


2 GMyg _ v2 
) erns? e 203 H (va — 2G Mns/TNs = Ves . (4.7) 


With Eqs. (4.4)-(4.6), the dark matter number density at rng can be obtained by 
integrating over all kinematically allowed velocities, v 2 \/2GMng/rng 


= alts) _ 


Na (rns) = = f 
Ma v>,/2GMyns/rNs 


For an homogeneous and isotropic (in the frame of the neutron star) axion phase space 
distribution far away from the neutron star and Vesce > Oy, one can find an approximate 
closed solution for na(rns), see Ref. [20]. We stress that their solution [and Eq. (4.7)] is only 
applicable if, far away from the neutron star, both the axion density and velocity distributions 
are both homogeneous and isotropic (in the frame of the neutron star) and if Vese > oy. A 
typical neutron star moves with speed ung of order oy relative to the Galactic rest frame, 
breaking the assumption of an isotropic velocity distribution in the frame of the neutron star. 

While the full six dimensional phase space contains all relevant information, for any 
realistic conversion location, \/2GMyg/rns >> {0v,Vesc}. Hence, almost the entirety of 
the axion’s speed comes from the infall. Thus, the distribution over v can approximately 
be treated as a delta-function 6(v — \/2GMyg/rng). In other words, we care about the 
distribution of incoming angles, rather than the distribution over incoming speeds. This 
allows us to define 


felens, /2GMns/rws, 9v, bv) 
Ov, n) = , 4.9 
Go oe) f a% felens, V2GMys/rys, 0, ¢',) (49) 


Thus armed, we can calculate the flux transfer for general axion distributions, including 
non-vanishing speeds of the neutron star relative to the galaxy. 


d°v ferns, v). (4.8) 


5 Comparison with 1D calculations in the Goldreich-Julian Model 


It is clear that for specific values 0 there can be significant variation between the 3D and 
1D versions of the axion-photon conversion. As the plasma is generated and shaped by 
the neutron star magnetic field, the changes in wp are driven by the changing B-field, and 
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so may also be correlated with 0. Here we explore what influence a more complete axion- 
photon conversion calculation has under realistic conditions. After a brief discussion of the 
Goldreich-Julian neutron star model in section 5.1, we will compare results for the axion- 
photon conversion in our full 3D calculation with results in the conventional 1D calculation 
in sections 5.2—5.3. 


5.1 The Goldreich-Julian Model 


To consider a realistic scenario, we turn to an analytic description of a neutron star. One 
commonly used model is the Goldreich-Julian (GJ) model [42]. The GJ model gives a simple 
analytic solution by requiring a charge density to satisfy the condition E-B = 0 at the surface 
of the neutron star in the presence of a strong, rotating dipole magnetic field given by 


3 
Brys = Bo (=) [cos Om cos ONS + sin Om sin Oys cos(Nt)] , 


B 3 
Boys = = (=) [cos 0m sin Ông — sin Om cos Ong cos(Nt)] , 


2 \rns 
Bo ( ro \* 
Bien = os (=) sin Ôm sin(Q¢) , (5.1) 


where Bo is the magnetic field strength at the surface of the neutron star (which is at located 
at ro), and Ôm is the misalignment angle between Byg and Q. The rotation vector of the 
neutron star is given by Q = NZng with Q = 27/T and T the neutron star spin period. We 
denote the polar angle of a given location with respect to the rotation axis zyg to be Ong. 
The charge density is given by 


22. - Bys 1 


; 5.2 
e 1 — 0?r? sin? Ong oe) 


nas(rys) = 


While we consider all charges to be electrons for simplicity, such that naj simply gives the 
number of electrons, our conversion calculation still applies to more complicated plasmas. 
The plasma frequency of an electron plasma is wp ~ \/47ane/mMe, giving 


An 1 
= B -2 . 
wp(rNs) y NS NS FT yan One’ (5.3) 
with 
a Bo / ro \° PO 
Bys : Zng = — | — [3 cos 6ys mh: fyg — cos Om] f (5.4) 
2 \rns 
Similarly to Ref. [14] we have also defined 
mM - yng = cos Ôm cos Ong + Sin Am sin Ong cos(Qt) . (5.5) 


It should be noted that it is unlikely that the GJ model holds over the entirety (or any) 
of the neutron star. The nature of neutron star atmospheres is highly uncertain, due to both 
the difficulty in modeling such large electromagnetic systems and the lack of observational 
data [43-46]. Asides from possible higher multipole components of the magnetic field [47], 
pair creation cascades at the poles may lead to significantly higher electron densities, as well 
as highly boosted electrons [45, 48, 49]. Our formalism allows for both y > 1 and ne £ naj, 
however, for specific examples we will remain within the GJ model. 
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Figure 2. Left panel: Flux transfer R as a function of incoming axion angles 6,, Ø» for the full 3D 
calculation. The neutron star has By = 10!4 G, 0m = 0.4 and T = 1s, and the ratio is calculated at 
the conversion zone located at Ong = 0.5, dng = 0.67 in the neutron star frame. We assume an axion 
mass of 25 eV with a coupling strength gay = 10714 GeV‘. Right panel: ‘Probability’ of axion to 
photon conversion ee as a function of incoming axion angles 6,, ¢, for a 1D calculation. The same 
neutron star and axion values are used. The black dashed line indicates the slice 6, = 7/3 that is the 
focus of Fig. 3. 


5.2 Directionality of emitted photons 


To see the impact of a 3D formalism, we can take the simplest example and ask: how does the 
flux transfer change for different incoming axion trajectories? As we are considering the GJ 
model, for all cases y = 1, meaning that üp = wp. In section 3 we derived the flux transfer R 
in a frame centered around the axion’s velocity (Z) and the direction of the B-field. However, 
as discussed above, the properties of the plasma are usually described in the frame of the 
neutron star, typically in spherical coordinates. To write down our X and y directions in 
terms of the incoming axion direction ĉ and magnetic field Ê we can use 


l a 
$ = B Z 5.6 

R= — Bus x Ê (5.6a) 
y =2xk (5.6b) 


For an example, we consider a neutron star with By = 10" G, 6,, = 0.4 and T = 1s. Taking 
some point on the sphere, Ong = 0.5, ns = 0.67, we can compare the flux transfer at the 
axion conversion surface as a function of the incoming axion angles, denoted by 0y, Ø» in the 
neutron star coordinate frame. Due to the dependence of the LO mode dispersion relation 
on @ and k, Eq. (2.19) indicates that the radius of conversion actually is slightly different for 
different incoming axion angles. However, the change in flux transfer due to small changes 
in conversion radius is subdominant compared to the impact of including aniosotropic effects 
of the medium. For a benchmark, we take the axion to have mass Mma = 25 eV and coupling 
Jay = 1071 GeV. Raising or lowering ma changes where the conversion zone is located: 
higher mass axions convert closer to the neutron star where magnetic fields are larger, usually 


=i = 


leading to greater flux transfers. If the axion mass is too large, however, it is greater than 
the plasma frequency everywhere outside the surface of the neutron star and no resonant 
conversion can occur. As the plasma frequency is set by the magnetic field strength of 
the neutron star, for each star, there will be a maximum m, which will allow for resonant 
axion-photon conversion. 

We plot the results of the 3D calculation, [the relativistic version of Eq. (3.3)], and of 
the 1D calculation, Eq. (3.4), in the right and left panels of Fig. 2, respectively. To avoid 
divergences when the WKB approximation breaks down, we use a simple cut-off estimated 
by the typical lengths scales over which the second order derivatives should come into play, 
derived in Appendix C 

2 p2 nr? le 

Reut-off = Jay BNS (Z) . (5.7) 

We see that the number of photons generated on a given trajectory can differ between the 

3D and the 1D formalism by several orders of magnitude. Furthermore, regions that were 

highly disfavoured in the 1D calculation can turn out to have maximal photon production 
when the aniosotropic nature of the medium is taken into account. 

To understand better the underlying physics behind the regions of enhanced and sup- 
pressed conversion, we can take a single slice of 6,,¢, space. In Fig. 3, we show a slice of 
constant 0y, in particular, 0, = 7/3. The top panel of Fig. 3 shows the ‘conversion probabil- 
ity’ for the 1D calculation and the flux transfer calculated with 3D effects, with the regulating 
cut-off flux transfer [Eq. (5.7)] shown by the gray dashed line. In the middle panel, we show 
the ratio R/ P Ra: up to three orders of magnitude enhancement can occur for specific axion 
trajectories. The bottom panel shows the relevant derivatives for each calculation, 0W,/0z 
in the 1D case and 0W,/0s for the 3D calculation. The order-of-magnitude discrepancies be- 
tween the 1D and 3D calculations of the flux transfer occur when 0W,/0z 4 OW,/Os. While 
there are still differences caused by correctly calculating the energy stored in the LO mode 
and the correct angular dependence, the primary difference between calculations comes from 
the difference between 8 and ĉ: when Ow,/Os ~ Ow,/Oz, the answers are similar up to the 
different angular dependence. The divergences in R, Pay (which are regulated by the cut- 
off) occur when Ow,/O(s,z) = 0, as indicated by the gray dashed line in the bottom panel 
of Fig. 3. In these examples, 6’ is negligible and so the cut-off is a stronger regulator. Note 
that the relationship between § and z does not depend on Ma: even for very non-relativistic 
axions which could convert far out from the NS there would still be significant differences in 
the flux transfer in a given direction. 

We can see that treating the LO mode as essentially a transverse mode (possibly with 
different dispersion, as in [23]) leads to a significant miscalculation of the axion-photon flux 
transfer. To have a reliable calculation of the axion-photon flux transfer for a given axion 
trajectory, one must take into account both the full energy stored in the LO mode, and, 
more importantly, the distinction between the direction of the axion’s momentum, ĉ, and the 
direction over which the energy density evolves, 8. 


5.3 Integrated power 


Asides from the directional dependence of the photons emitted from a given point, one 
may also be concerned with the total power emitted over some portion of the neutron star. 
Studies of ray-tracing suggest that the photons bend relatively quickly as plasma frequency 
drops [23, 24]. The momentum acquired by the photon as it transitions to vacuum must 
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Figure 3. Top panel: Axion-photon flux transfer R as a function of incoming axion azimuthal angle 
v for the full 3D calculation (green) and 1D calculation (blue). For the 1D calculation we use the 
‘conversion probability’ instead of the flux transfer. Middle panel: Ratio of the full 3D axion-flux 
transfer R to the 1D conversion probability. The flux is enhanced by up to three orders of magnitude 
with resepect to the 1D calcuation. Bottom panel: Rate of change of the plasma frequency wp 
with respect to s (green) and z (blue). Axion photon conversion diverges when Ow,/O(z,s) = 0, the 
gray dashed line. The neutron star has By = 10/4 G, 0m = 0.4 and T = 1s, and the quantities are 
calculated at the conversion zone located at Ons = 0.5,¢@ng = 0.67 in the neutron star frame. We 
assume an axion mass of 25 ueV with a coupling strength gay = 10714 GeV! and incoming polar 


angle 0, = 7/3. 


come from the gradient of the medium, so photons that begin mildly relativistic would 
become approximately parallel to the change in Ww». Because of this, one can make a rough 
estimate that the total power exiting the neutron star will end up on the same trajectory, 
determined by Væp. We stress however that the axion can be reasonably relativistic, and so 
one should do a full ray-tracing calculation for reliable results. 
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Figure 4. Axion photon conversion integrated over axion phase space, Ry, as a function of incoming 
neutron star polar angle Oys for 3D (green) and 1D (blue) calculations. The 3D calculation uses the 
full flux transfer of axions to photons, whereas the 1D case only considers the ‘conversion probability’ 
interpreted from a ‘Schrédinger-like’ equation. We choose a neutron star with Bp = 10 G, Om = 0.4 
and T = 1s. The functions are calculated at the conversion zone located at an azimuthal angle 
ons = 7 in the neutron star frame. We assume an axion mass of 25 weV with a coupling strength 
Jay = 10714 GeV~! as well as an isotropic velocity distribution of the axions. We mark the region 
where the conversion radius would be less than the neutron star radius rọ with the red dashed lines 
labeled ‘Inside NS’. 


To calculate the power radiated from a given point on the neutron star rys, we must 
integrate over the incoming axion angles 6», Ø», weighted with an appropriate phase space 
distribution. To focus on the difference caused by properly treating the aniostropy of the 
medium, we will neglect the small differences in conversion radii for different axion trajectories 
and use a single radius of conversion given by the transverse photon dispersion, Ma = Wp. 
The variation due to choosing different radii seems to be at most a factor of 2 when one 
considers the longitudinal photon dispersion (w = wp) instead, and does not change the 
general relationship between the two calculations. Similarly, we neglect contributions from 
6’, which again gives a subdominant contribution for the cases we consider. The distribution- 
weighted axion-photon flux transfer Ry is given by 


Rw (rns) = J e%f5lens, 80. o) RENS: bor o), (5.8) 


where fs is the five dimensional phase space defined in Eq. (4.9). For the 1D case we use the 
1D ‘conversion probability’, Eq. (3.4) instead of R. 

We can now show Ry for a slice along the neutron star conversion surface assuming an 
isotropic axion distribution. We keep dng fixed at 7 and vary Ong to derive Fig. 4. The struc- 
ture generated in a slice through Ong is in general behaviour independent of the specifically 
chosen @ng. Integrating over all photon directions generally smooths out differences between 
the 3D and 1D calculations. However, there are still noticeable differences. Overall, the 3D 
calculation leads to more energy stored in the EM fields, leading to more photons radiated 
from the neutron star. As the 1D calculation neglected energy stored in the non-transverse 
components of the Hamiltonian, we would naturally expect to see additional ‘conversion’ just 
from correctly counting the number of photons generated. This enhancement is particularly 


a 


noticeable around the neutron star ‘throat’, i.e. the region where the conversion radius van- 
ishes as Q - Bys — 0, just outside the red boarders in Fig. 4. Inside the neutron star itself, 
the plasma density no longer follows the GJ model; given the extreme densities inside a neu- 
tron star the plasma frequency is much larger than the axion masses of interest everywhere 
inside a neutron star. Thus, no propagating photons can be produced from axion-photon 
conversion inside a neutron star. One would expect that issues of photon curvature, as ex- 
plored in Ref. [23], as well any concerns about gravity bending the axion path would be most 
pronounced in the throats, making this region generally unreliable. One would likely need 
to do numerical calculations for a robust solution including a gravitationally bending axion 
trajectory. 

Having a similar total conversion rate for isotropic dark matter infall is not surprising, 
as ultimately the length scales and magnetic fields that give the general order of conversion 
are the same. However, where the axions are converted, and in what directions the pho- 
tons emerge, is significantly impacted. The plasma effectively lenses the outgoing photons, 
which when combined with the limited line of sight and Doppler shifting of light can give a 
complicated, time dependent signal. Thus any aspect which changes the location and direc- 
tion of conversion will impact axion searches significantly. Because of this we focus on the 
changes to the power and directionality of a given point on the neutron star, as the total 
power radiated from a neutron star is not usually relevant to an individual observation of an 
individual neutron star. In addition, for compact clumps of axions, such as axion stars [50- 
55] or miniclusters [56-61], one would expect a much stronger dependence in the resulting 
signal depending at what angle the compact object hits the neutron star. Thus, changes 
to the production location and direction can have a disproportionate effect on observational 
constraints. 


6 Conclusions 


The axion is one of the premier candidates for dark matter. It displays unique wavelike 
phenomenology requiring both novel techniques in the lab and astrophysical searches. Reso- 
nant conversion to photons in the atmospheres of neutron stars is one such possible method 
for discovering the axion. Neutron stars are complicated systems and therefore calculating 
the projected signal is extremely challenging. Here we have focused on one piece of the 
larger puzzle, and together with recent advances in ray-tracing help lay the road to a robust 
prediction for what might be observed on Earth. 

To this end, we have calculated the production of photons via axions in strongly mag- 
netised, anisotropic environments. The prototypical example of such an environment is a 
neutron star, but our calculation is generic. Previous calculations have not explored the 
anisotropic nature of the medium and have focused only on the transverse part of pho- 
ton modes. We show, however, that the mixing of transverse and longitudinal modes into 
Langmuir-O modes causes the direction over which the amplitude of the electromagnetic 
fields evolve to be different from the direction of the axion. We find that rather than evolv- 
ing in the direction of propagation, the waves evolve also in the transverse direction. This 
directional difference can lead to orders of magnitude change in the flux transfer R in a given 
direction. 

Furthermore, taking only the transverse component of the stored energy into account 
neglects a significant part of the stored energy, or in other words, miscounts the number 
of photons that are converted. Taking only the transverse component essentially confuses 
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a classical, Schrédinger-like solution with the quantum mechanical conversion probability. 
We show that taking both features into account significantly modifies the probability of 
conversion over a given axion-trajectory. While a truly complete solution would require 
numerically solving the equations of motion in a varying medium, we provide simple analytical 
formulae which can be used to estimate the radio signals generated in neutron stars by axions. 
This will allow the robustness of axion limits to be increased in the case of non-observation, 
such as Refs. [25-28], or allow one to estimate the neutron star and axion properties in the 
case of a discovery. 
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A Longitudinal component 


In the main text, we solve the axion-Maxwell equations by eliminating the longitudinal terms 
from the equations. However, this does not mean that the longitudinal modes do not exist. 
To solve for the longitudinal component excited by the axion, we can rewrite Eq. (2.8) to 
find i 
aoa — iT cos 0 sin 6, — Ww gaya Bys cos 0 
7 : 


E, = (A.1) 


w2 — 52 cos2 0 


First we should consider the case without an axion to find the normal electromagnetic modes 
of the system. In the absence of the axion, for a homogeneous medium one finds that 

2 
aoe — w? cos 0 sin 0 


p 
w2 — we cos? 0 By, ue) 


E, = 


allowing us to define the propagating Langmuir-O (LO) mode, which becomes an ordinary 
transverse wave when 0 = 7/2 and longitudinal when 6 = 0 in the non-relativistic limit [33], 


ž —1 _  @cosdsind | 
Eo = k 2 ) , (A.3) 


= Z 
1 ws cos? 0 sin? 6 We — w2 cos2 0 
+ (2—02 cos2 0)? 


where we have inserted a factor of —1 so that ÊLo = Byg in the non-relativistic limit. 
Including the axion current, and knowing that we will solve for Ey directly, we can neglect 
the subdominant derivative terms to find 


-2 . 2 
A w-cos@sind -~ w° cos 0 
E, ~ z E, Jay BNså . A.4 
Á w2? — w cos20 ” w2— w2 cos? 0 a (AA) 
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The first term is just the longitudinal component of the LO mode, with the second cor- 
responding to direct axion mixing. In the limit of a perfectly aligned magnetic field with 
0 = {0,7} Eq. (A.4) reduces to the expression for a pure longitudinal excitation in a lossless 
plasma [31, 62], 


= B 
Ë, Jay Ne (A.5) 
1-3 


We can see that the two different terms in FE, have very different propagation behaviour. 
After some distance (given by the conversion length L) Ey, will be out of phase with a. 
Notice that the terms proportional to à are not generally propagating: as soon as the axion 
mixing goes away, so to does this contribution to Ey. 

Thus, for the purpose of solving for an outgoing radio signal from the neutron star, the 
mode that the axion couples to is the LO mode, which will adiabatically transform into a 
transverse mode as it exits the system [33]. 


B Stationary phase approximation 


To provide an analytic expression for axion-photon conversion, in Eq. (2.15) in the main text 
we derived a Schroédinger-like equation 


OE 1 “a x lwe | 
Fg = Dp (Ma — Bp — IRD) By a Bws B) 
which has solution 
: S 1 w2 Pa n = L Em 
m=- fds E oäByse Si Deika klm), (B2) 


Rather than a numeric solution, we can use the stationary phase approximation to find 
an analytic solution. The premise of the stationary phase approximation is that highly 
oscillatory integrands tend to cancel when integrated over many oscillation periods, and 
that the dominant contribution therefore often comes from stationary points, say so, where 
the phase, f(s), is locally constant, i.e. f'(sọ) = 0. In the standard form of the stationary 
phase approximation, the non-oscillatory factor of the integrand, say g(s), is approximated as 
constant around so, and the phase is expanded to quadratic order: f(s) ~ f(so)+4f”(so)(s— 
so)’, and the limits of integration are taken to too. For s > so, this gives the approximation, 


et if(s if” leo) (g s if(s sign| f” (so)lir 2T 
[eats if(s ~ g(so)e if ( Fa ds e 2 ™ (5-50)? = g( soJe! o)+sign[f” (so)lin/4 F 


(B.3) 
The quality of the stationary phase approximation of course depends on the accuracy 
of the assumptions that it relies on. It’s convenient to define a (squared) ‘conversion length’ 


as 
T 


|£” (so) 


In the limit where ,/|f”(so)| is larger than all other scales in the problem, the conversion 
length is small, and the stationary phase approximation is expected to be very accurate. 
However, in many physical situations there are no large hierarchies of scales, and one needs 
to consider corrections to the standard stationary phase formula. Moreover, extending the 


i= (B.4) 
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integration range to +00 is not always justified, and the full integral that one wants to 
evaluate may only be well-defined in a small region around the stationary point, such as for 
our Eq. (B.2). 

In the context of relativistic axion-photon conversion, these issues were recently dis- 
cussed in Ref. [63], where corrections to Eq. (B.3) were derived from considering a finite in- 
tegration range, and accounting for contributions neglected in the standard stationary phase 
approximation. For the purpose of our analysis, the main point is that when one considers 
the integration range [so — As, so + As], the stationary phase approximation gives 


s ; : sot+tAs ifls 
| ds g(s)eF) = g(sp)e'F 0) l ds e F (6-30)? (B.5) 
0 so—As 
= g(so)etf(s0)+sianLF"(so)]ia/4./9 7 Be ( mA) , (B.6) 


where Err denotes the error function. The norm of the error function grows linearly with 
as for small arguments, and performs slowly damped oscillations around unity for as > 1. 
This in particular implies that the standard formula only applies when it is justified to take 
As > L. Equation (B.6) applies for any ratio of As/L, but simplifies considerably when 
Az < L, and the integral becomes 


| f ds g(s) ®| = 2g(s0)|Az, (B.7) 


up to corrections to cubic order in Az/L. This implies that the formula (B.3) can be regulated 
and modified in a very simple way, by replacing L > /2Az, when it is justified to consider 
Az < L, for example when the WKB approximation breaks down, or the assumption of 
constant 0 fails. For a more detailed discussion of how to evaluate also the contributions 
from |z| > Az, see Ref. [63]. 

The most important factor comes from f”(sq). While we have so far considered 6 to be 
constant, if it is slowly changing the largest impact will be felt in the conversion length [23] 
as the impact of dephasing the axion and photon is much greater than the small changes to 
the definition of y. We can then see that 


CE Me ee ot Gp, sin? 0 a? (w? — w?)6' sin 0 cos 0 
S- ( og [m2 - £03] ) = p. s - a (B.8) 
k (1 — z4 cos? 6) wk (1 — 5 cos? 6) 


We now apply the standard stationary phase approximation of Eq. (B.3) to the electric 
field of Eq. (B.2), which gives 


= ol wog 520! 
2k pin, + we tan 0” p 


£ T 50 qe! - 
Ey[so + L/2] = n elo” ds D/2 JayăBys (B.9) 


where we have neglected an overall phase and where œ, = OW,/0s. This is the analytical 
form of axion-induced E, that we will use in our subsequent analysis. The stationary point 


condition (cf. f’(s9) = 0) translates into the ‘resonance condition’: 


2, 2 
maw 


p(s)? = (B.10) 


m2 cos? 6 + w? sin? 6 ` 


= 9] = 


In the non-relativistic limit, w — mq, and the resonance condition simply becomes p = Ma. 
The conversion length evaluates to 


in 0 k 
= aa (B.11) 
E Dol + wep w26! 
Pp w2 tan@~ Pp 


For most axion trajectories, this conversion length is smaller than the length scales associated 
with the gravitational curvature of the axion trajectory, the curvature of the photon path 
due to the anisotropic medium, and the scale over which the magnetic field changes. When 
this is the case, the standard stationary phase approximation applies. However, this is not 
always the case, and in some cases the conversion length can get significantly enhanced, as 
we now discuss. 


C Regulating divergences 


The flux transfer R, given by Eq. (3.3), was calculated assuming the WKB approximation 
holds: i.e., that the axion momentum is much larger than the first derivatives of By, which are 
larger than the second derivatives. However, R diverges in the limit that the first derivative 
vanishes (i.e., the conversion length becomes infinite). Such a divergence must be regulated 
for two reasons. The first is the breakdown of the WKB approximation, and the second is 
the breakdown of the assumption that the conversion happens over length scales smaller than 
other relevant ones in the neutron star, for example the radius of photon bending due to the 
changing media. 

In the first case, it is evidently untrue that the vanishing first derivatives are larger 
than the second derivatives. Rather than attempting to solve the full 3D equations by hand 
with higher derivatives, as Ez, Ey and E, all appear with their second derivatives, we can 
use a simple estimate to institute a cut off on the diverging flux transfers. In the 1D case, it 
was shown in Ref. [14] that an order-of-magnitude estimate could be obtained by using the 
simple formula 


, (C.1) 


where L is the conversion length. The assumption of coherent conversion requires that the 
change in phase between two locations is small. Thus, we can estimate L by checking the 
distance over which k changes. We can approximate k changing over some length L by 
writing 

k(L) = ko + kL + kaL’, (C.2) 
where destructive interference occurs when ôk = k— ko satisfies dkL ~ m (as the phase enters 
via e(kL—“t) As we are concerned with the case where kıL S koL?, we can neglect the linear 


term and see that 
L? ~ r/k. (C.3) 


As the only change in k comes from the changing plasma frequency W,, in the limit of 
vanishing first derivatives kk2 = @p)W, 2, where similarly 


Wp(L) ~ Gpo + Gp2L. (C.4) 


To get a conservative estimate for L, we can use that at second order, all directions may 
play arole. Thus while the second derivatives with respect to some directions might be small 


= 9) = 


or vanishing, we know that the radial direction will have a power law dependence. Thus, 
assuming that the r direction is the relevant one will give a fairly conservative estimate for 
the total conversion probability. We stress that this is not a full solution of the 3D Maxwell 
equations, but rather a conservative way of regulating our WKB formalism. Taking the 
power law dependence of p, we then find that 


L= o : (C.5) 
Wp 

Putting in some typical numbers, we see that the typical conversion length should increase by 

about an order of magnitude, when compared to purely radial estimates of Ref. [14], leading 

to a hundredfold increase in the conversion probability. 

The second case relates to the overall trajectory of photons in the medium, for example, 
that the photon bends and so dephases with the axion, or other curvature effects of the 
medium [23]. While such a constraint may provide a more stringent limit, particularly in 
the neutron star throat, to calculate such effects in full requires ray tracing calculations [20, 
23, 24], which are beyond the scope of this work. When these effects are considered, one 
should use the stricter of the two cut-offs. As we are only interested in demonstrating the 
influence of a 3D, aniosotropic calculation of the flux transfer rather than a start to finish 
signal calculation for observational purposes, we can simplify matters and just use a cut-off 
based on the second derivatives of the E-field, 


rr 2/3 
Reut-off = Ion Bis (2 5) . (C.6) 
c 


References 
1| R. D. Peccei and H. R. Quinn, CP Conservation in the Presence of Instantons, Phys. Rev. 
Lett. 38 (1977) 1440. 


2| R. D. Peccei and H. R. Quinn, Constraints Imposed by CP Conservation in the Presence of 
Instantons, Phys. Rev. D 16 (1977) 1791. 


3] S. Weinberg, A New Light Boson?, Phys. Rev. Lett. 40 (1978) 223. 


4| F. Wilczek, Problem of Strong P and T Invariance in the Presence of Instantons, Phys. Rev. 
Lett. 40 (1978) 279. 


5| J. Preskill, M. B. Wise and F. Wilczek, Cosmology of the Invisible Axion, Phys. Lett. B 120 
(1983) 127. 


6] L. Abbott and P. Sikivie, A Cosmological Bound on the Invisible Axion, Phys. Lett. B 120 
(1983) 133. 


7] M. Dine and W. Fischler, The Not So Harmless Axion, Phys. Lett. B 120 (1983) 137. 


8] L. D. Duffy and K. van Bibber, Azions as Dark Matter Particles, New J. Phys. 11 (2009) 
105008 [0904.3346]. 


9] P. Arias, D. Cadamuro, M. Goodsell, J. Jaeckel, J. Redondo and A. Ringwald, WISPy Cold 
Dark Matter, JCAP 06 (2012) 013 [1201.5902]. 


[10] I. G. Irastorza and J. Redondo, New Experimental Approaches in the Search for Axion-Like 
Particles, Prog. Part. Nucl. Phys. 102 (2018) 89 [1801 .08127]. 


[11] Y. K. Semertzidis and S. Youn, Axion Dark Matter: How to Detect it?, 2104.14831. 


— 93 — 


12 


13 


14 


15 


16 


17 


18 


19 


20 


21 


22 


23 


24 


25 


26 


27 


28 


29 


30 


31 


F. Chadha-Day, J. Ellis and D. J. E. Marsh, Axion Dark Matter: What is it and Why Now?, 
2105.01406. 


M. S. Pshirkov and S. B. Popov, Conversion of Dark Matter Axions to Photons in 
Magnetospheres of Neutron Stars, J. Exp. Theor. Phys. 108 (2009) 384 [0711.1264]. 


A. Hook, Y. Kahn, B. R. Safdi and Z. Sun, Radio Signals from Axion Dark Matter Conversion 
in Neutron Star Magnetospheres, Phys. Rev. Lett. 121 (2018) 241102 [1804.03145]. 


D. Lai and J. Heyl, Probing Azions with Radiation from Magnetic Stars, Phys. Rev. D 74 
(2006) 123003 [astro-ph/0609775]. 


F. P. Huang, K. Kadota, T. Sekiguchi and H. Tashiro, Radio Telescope Search for the Resonant 
Conversion of Cold Dark Matter Azions from the Magnetized Astrophysical Sources, Phys. Rev. 
D 97 (2018) 123001 [1803 .08230]. 


B. R. Safdi, Z. Sun and A. Y. Chen, Detecting Axion Dark Matter with Radio Lines from 
Neutron Star Populations, Phys. Rev. D 99 (2019) 123021 [1811.01020]. 


T. D. P. Edwards, M. Chianese, B. J. Kavanagh, S. M. Nissanke and C. Weniger, Unique 
Multimessenger Signal of QCD Axion Dark Matter, Phys. Rev. Lett. 124 (2020) 161101 
[1905 .04686]. 


R. A. Battye, B. Garbrecht, J. I. McDonald, F. Pace and S. Srinivasan, Dark Matter Axion 
Detection in the Radio/mm- Waveband, Phys. Rev. D 102 (2020) 023504 [1910.11907]. 


M. Leroy, M. Chianese, T. D. P. Edwards and C. Weniger, Radio Signal of Axion-Photon 
Conversion in Neutron Stars: A Ray Tracing Analysis, Phys. Rev. D 101 (2020) 123003 
[1912.08815]. 


J. H. Buckley, P. S. B. Dev, F. Ferrer and F. P. Huang, Fast Radio Bursts from Axion Stars 
Moving Through Pulsar Magnetospheres, Phys. Rev. D 103 (2021) 043015 [2004.06486]. 


S. Nurmi, E. D. Schiappacasse and T. T. Yanagida, Radio Signatures From Encounters Between 
Neutron Stars and QCD-Azion Minihalos Around Primordial Black Holes, 2102.05680. 


S. J. Witte, D. Noordhuis, T. D. P. Edwards and C. Weniger, Azion-Photon Conversion in 
Neutron Star Magnetospheres: The Role of the Plasma in the Goldreich-Julian Model, 
2104.07670. 


R. A. Battye, B. Garbrecht, J. I. Mcdonald and S. Srinivasan, Radio Line Properties of Axion 
Dark Matter Conversion in Neutron Stars, 2104.08290. 


J. W. Foster, Y. Kahn, O. Macias, Z. Sun, R. P. Eatough, V. I. Kondratiev, W. M. Peters, 

C. Weniger and B. R. Safdi, Green Bank and Effelsberg Radio Telescope Searches for Axion 
Dark Matter Conversion in Neutron Star Magnetospheres, Phys. Rev. Lett. 125 (2020) 171301 
[2004. 00011]. 


J. Darling, New Limits on Axionic Dark Matter from the Magnetar PSR J1745-2900, 
Astrophys. J. Lett. 900 (2020) L28 [2008.11188]. 


J. Darling, Search for Axionic Dark Matter Using the Magnetar PSR J1745-2900, Phys. Rev. 
Lett. 125 (2020) 121103 [2008 . 01877]. 


R. A. Battye, J. Darling, J. McDonald and S. Srinivasan, Towards Robust Constraints on 
Axion Dark Matter Using PSR J1745-2900, 2107 .01225. 


G. Raffelt and L. Stodolsky, Mixing of the Photon with Low Mass Particles, Phys. Rev. D 37 
(1988) 1237. 


P. Sikivie, Experimental Tests of the Invisible Axion, Phys. Rev. Lett. 51 (1983) 1415. 
(Erratum: Phys.Rev.Lett. 52, 695 (1984)]. 


A. J. Millar, J. Redondo and F. D. Steffen, Dielectric Haloscopes: Sensitivity to the Axion 
Dark Matter Velocity, JCAP 10 (2017) 006 [1707.04266]. [Erratum: JCAP 05, E02 (2018). 


— 94 — 


32 


33 


34 


35 
36 


37 


38 


39 


40 


41 


42 
43 


44 


45 


46 


47 


48 


49 


50 
51 


52 
53 


A. V. Gurevich, V. S. Beskin and Y. N. Istomin, Physics of the Pulsar Magnetosphere. 
Cambridge University Press, 1993, 10.1017/CBO9780511564673. 


M. Gedalin, D. B. Melrose and E. Gruman, Long Waves in a Relativistic Pair Plasma in a 
Strong Magnetic Field, Phys. Rev. E 57 (1998) 3399. 


A. Caputo, A. J. Millar and E. Vitagliano, Revisiting Longitudinal Plasmon-Azion Conversion 
in External Magnetic Fields, Phys. Rev. D 101 (2020) 123004 [2005 . 00078]. 


A. Bers, Space-Time Evolution of Plasma Instabilities - Absolute and Convective, 1983. 


A. Bers, Note on Group Velocity and Energy Propagation, American Journal of Physics 68 
(2000) 482. 


G. Raffelt, G. Sigl and L. Stodolsky, Quantum Statistics in Particle Mixing Phenomena, Phys. 
Rev. D 45 (1992) 1782. 


M. S. Alenazi and P. Gondolo, Phase-space Distribution of Unbound Dark Matter Near the 
Sun, Phys. Rev. D 74 (2006) 083518 [astro-ph/0608390]. 


A. K. Drukier, K. Freese and D. N. Spergel, Detecting Cold Dark Matter Candidates, Phys. 
Rev. D 33 (1986) 3495. 


J. D. Lewin and P. F. Smith, Review of Mathematics, Numerical Factors, and Corrections for 
Dark Matter Experiments Based on Elastic Nuclear Recoil, Astropart. Phys. 6 (1996) 87. 


K. Freese, M. Lisanti and C. Savage, Colloquium: Annual Modulation of Dark Matter, Rev. 
Mod. Phys. 85 (2013) 1561 [1209.3339]. 


P. Goldreich and W. H. Julian, Pulsar Electrodynamics, Astrophys. J. 157 (1969) 869. 


C. Kalapotharakos, D. Kazanas, A. Harding and I. Contopoulos, Toward a Realistic Pulsar 
Magnetosphere, Astrophys. J. 749 (2012) 2 [1108.2138]. 


B. Cerutti and A. Beloborodov, Electrodynamics of Pulsar Magnetospheres, Space Sci. Rev. 
207 (2017) 111 [1611.04334]. 


A. A. Philippov and A. Spitkovsky, Ab-Initio Pulsar Magnetosphere: Particle acceleration in 
Oblique Rotators and High-energy Emission Modeling, Astrophys. J. 855 (2018) 94 
[1707 P 04323]. 


M. G. Baring, Z. Wadiasingh, P. L. Gonthier, A. K. Harding and K. Hu, The Mysterious 
Magnetospheres of Magnetars, PoS HEASA2019 (2021) 036 [2012.10815]. 


A. V. Bilous, A. L. Watts, A. K. Harding, T. E. Riley, Z. Arzoumanian, S. Bogdanov, K. C. 
Gendreau, P. S. Ray, S. Guillot, W. C. G. Ho and D. Chakrabarty, A NICER View of PSR 
J0030+0451: Evidence for a Global-scale Multipolar Magnetic Field, The Astrophysical Journal 
887 (2019) L23. 


T. Erber, High-Energy Electromagnetic Conversion Processes in Intense Magnetic Fields, Rev. 
Mod. Phys. 38 (1966) 626. 


A. K. Harding and D. Lai, Physics of Strongly Magnetized Neutron Stars, Reports on Progress 
in Physics 69 (2006) 2631. 


D. J. Kaup, Klein-Gordon Geon, Phys. Rev. 172 (1968) 1331. 


R. Ruffini and S. Bonazzola, Systems of Selfgravitating Particles in General Relativity and the 
Concept of an Equation of State, Phys. Rev. 187 (1969) 1767. 


I. I. Tkachev, On the Possibility of Bose Star Formation, Phys. Lett. B 261 (1991) 289. 


P.-H. Chavanis, Mass-Radius Relation of Newtonian Self-Gravitating Bose-Einstein 
Condensates with Short-Range Interactions: I. Analytical Results, Phys. Rev. D 84 (2011) 
043531 [1103.2050]. 


— 95 — 


54 


55 


56 
57 


58 


59 


60 


61 


62 


63 


L. Visinelli, S. Baum, J. Redondo, K. Freese and F. Wilczek, Dilute and Dense Axion Stars, 
Phys. Lett. B 777 (2018) 64 [1710.08910]. 


A. Prabhu and N. M. Rapidis, Resonant Conversion of Dark Matter Oscillons in Pulsar 
Magnetospheres, JCAP 10 (2020) 054 [2005 .03700]. 


C. J. Hogan and M. J. Rees, Axion Miniclusters, Phys. Lett. B 205 (1988) 228. 


E. W. Kolb and I. I. Tkachev, Axion Miniclusters and Bose Stars, Phys. Rev. Lett. 71 (1993) 
3051 [hep-ph/9303313]. 


E. W. Kolb and I. I. Tkachev, Nonlinear Axion Dynamics and Formation of Cosmological 
Pseudosolitons, Phys. Rev. D 49 (1994) 5040 [astro-ph/9311037]. 


E. W. Kolb and I. I. Tkachev, Large Amplitude Isothermal Fluctuations and High Density Dark 
Matter Clumps, Phys. Rev. D 50 (1994) 769 [astro-ph/9403011]. 


B. Eggemeier, J. Redondo, K. Dolag, J. C. Niemeyer and A. Vaquero, First Simulations of 
Axion Minicluster Halos, Phys. Rev. Lett. 125 (2020) 041301 [1911.09417]. 


T. D. P. Edwards, B. J. Kavanagh, L. Visinelli and C. Weniger, Transient Radio Signatures 
from Neutron Star Encounters with QCD Axion Miniclusters, 2011.05378. 


M. Lawson, A. J. Millar, M. Pancaldi, E. Vitagliano and F. Wilczek, Tunable Axion Plasma 
Haloscopes, Phys. Rev. Lett. 123 (2019) 141802 [1904.11872]. 


M. C. D. Marsh, J. H. Matthews, C. Reynolds and P. Carenza, The Fourier Formalism for 
Relativistic Axion-Photon Conversion, with Astrophysical Applications, 2107 . 08040. 


— 26 — 


